#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>

#define max_iteration 1E3

double f( double x ){
	return x*x - 1;
}

/*
 * findroot: 在[a,b]区间内寻找 func(x)=0　的根，并返回根的值
 * 算法：二分法
 */
double bisection(double (*func)(double), double a, double b, double precision){

	ofstream fout( "bisec.gnu" );
        fout<< " set term gif animate \n set output \"bisec.gif\" \n ";
        fout<< " set samples 10000 \n ";
        fout<< " plot [0:"<<b<<"] x*x-1 w l t \"y=x^2\", 0 w l lc \"black\" notitle "<<endl;

        double ya=func(a), yb=func(b);
        double c=(a+b)/2, yc;
        int i;

        if( ya* yb > 0){//如果 func(a) * func(b) >0, 则[a,b]区间内不一定有解，拒绝任务
                cout<<"error: invalid input [a,b] for findroot: f(a)*f(b)>0 "<<endl;
                exit(1);
        }

        for(i=0;i<max_iteration;i++){

                c = (a+b)/2;//区间中点
                yc = func(c);
		for(int j=0;j<10;j++){
			fout<< " set arrow "<<i+1<<" nohead from "<<c<<",0 to "<<c<<","<<yc<<" nohead lw 2 lc \"red\" "<<endl;
			fout<< " set arrow "<<10*(i+1)<<" nohead from "<<a<<",0 to "<<b<<",0 nohead lw 2 lc \"red\" "<<endl;
			fout<< " set title \" Point "<<i+1<< "\" font \"Arial,20\" "<<endl;
			fout<< " replot "<<endl;
			fout<< " unset arrow "<<i+1<<"\n";
		       	fout<< " unset arrow "<<10*(i+1)<<endl;
			fout<<" replot \n";
		}

                if( yc * ya <= 0 ){//如果 yc*ya<=0，说明[a,c]中有一个根
                        b=c;
                }
                else{//如果 yc*ya>0，说明[c,b]中有一个根
                        a=c;
                }
                if( fabs(a-b) < precision ){
                        cout<<"\n iteration rounds = "<<i+1<<endl;
			fout.close();
                        return (a+b)/2;
                }
        }
        cout<<"bisection method: failed to find a root, after "<<max_iteration<<" iterations.\n";
        exit(1);
}

void fdf( double x, double & y, double & dy ){
	y = x*x - 1;
	dy = 2*x;
};

/*
 * Newton's method: given a function (*func), which gives value and also the 1st derivative,
 * a suspicious zone [a,b], it finds a root in the zone within precision, and returns the root.
 */
double Newton_findroot(void (*func)(double, double &, double &), double x0, double precision){

        double x, dx, y, dy;
        int i;

	ofstream fout("newton.gnu");
	fout<< " set term gif animate \n set output \"newton.gif\" \n ";
	fout<< " set samples 10000 \n ";
	fout<< " plot [0:"<<x0<<"] x*x-1 w l t \"y=x^2\", 0 w l lc \"black\" notitle "<<endl;

        x = x0;

        for(i=0;i<max_iteration;i++){
                func(x, y, dy);
		fout<< " set arrow "<<i+1<<" nohead from "<<x<<",0 to "<<x<<","<<y<<" nohead lw 2 lc \"red\" "<<endl;
		fout<< " set title \" Point "<<i+1<< "\" font \"Arial,20\" "<<endl;
                if( fabs(dy)<1E-9 ){
                        cout<<" error: f'=0 in Newton's root-finding method. \n";
                        exit(1);
                }
                dx = y/dy;
		double tempx = x;
                x -= dx;
	
		for(int j=0;j<10;j++){
			fout<< " set arrow "<<10*(i+1)<<" from "<<tempx<<","<<y<<" to "<<x<<","<<0<<" nohead lw 2 lc \"red\""<<endl;
			fout<< " replot "<<endl;
			fout<< " unset arrow "<<10*(i+1)<<" \n replot \n ";
		}
		fout<< " set arrow "<<10*(i+1)<<" from "<<tempx<<","<<y<<" to "<<x<<","<<0<<" nohead lw 2 lc \"red\""<<endl;
	
                if( fabs(dx) < precision ){
                        cout<<"\n iteration rounds = "<<i+1<<endl;
			fout.close();
                        return x;
                }
        }
        cout<<" Newton's method: after "<<max_iteration<<" steps, failed to find a root.\n";
        exit(1);
}

int main(){

	double a = 0, b = 4, precision = 1E-9;

	double x_bisec = bisection( f, a, b, precision );

	double x_newton = Newton_findroot( fdf, b, precision);

	return 0;
}
